
##load data
x.vars = c('non.white','hispanic','college.grad','male')
money.dat = read.csv('money_dat.csv')
attribution.dat = read.csv('attribution_dat.csv')
social.dat = read.csv('social_dat.csv')
physical.dat = read.csv('physical_dat.csv')

##limit data to subjects in the first row
money.dat = money.dat[money.dat$disbelief == 0,]
attribution.dat = attribution.dat[attribution.dat$disbelief == 0,]
social.dat = social.dat[social.dat$disbelief == 0,]
physical.dat = physical.dat[physical.dat$disbelief == 0,]

##if limiting, it can screw up blocks because small sessions might have nobody randomly assigned to first row
##so need to only use only blocks with treatment and control
block.assignments = table(money.dat$Block.ID,money.dat$treatment.assignment)
block.assignments = block.assignments[block.assignments[,1] !=0 & block.assignments[,2] !=0,]
complete.blocks = as.numeric(row.names(block.assignments))
money.dat = money.dat[money.dat$Block.ID%in%complete.blocks,]

block.assignments = table(attribution.dat$Block.ID,attribution.dat$treatment.assignment)
block.assignments = block.assignments[block.assignments[,1] !=0 & block.assignments[,2] !=0,]
complete.blocks = as.numeric(row.names(block.assignments))
attribution.dat = attribution.dat[attribution.dat$Block.ID%in%complete.blocks,]

block.assignments = table(social.dat$Block.ID,social.dat$treatment.assignment)
block.assignments = block.assignments[block.assignments[,1] !=0 & block.assignments[,2] !=0,]
complete.blocks = as.numeric(row.names(block.assignments))
social.dat = social.dat[social.dat$Block.ID%in%complete.blocks,]

block.assignments = table(physical.dat$Block.ID,physical.dat$treatment.assignment)
block.assignments = block.assignments[block.assignments[,1] !=0 & block.assignments[,2] !=0,]
complete.blocks = as.numeric(row.names(block.assignments))
physical.dat = physical.dat[physical.dat$Block.ID%in%complete.blocks,]


cat('\n \n Starting Experiment 2 Estimation, Non-Skeptical Subjects Only \n \n') 

##outmatrix
outmat = matrix(nrow = 4, ncol = 6)
rownames(outmat) = c('allocation','attribution','social.perceptions','physical.perceptions')
colnames(outmat) = c('int.mean','seg.mean','ate','p.value','cohens.d','N')

###allocation analysis
##RI analysis
money.Xs = as.matrix(money.dat[,x.vars])  ##get necessary x variables
##generate RI permutations and probabilities
money.ri.perms = genperms(money.dat$treatment.assignment,
                        blockvar = money.dat$Block.ID)
money.ri.probs = genprob(money.ri.perms)

##estimate treatment effect and inference
ate = estate(Y = money.dat$allocated.diff,
             Z= money.dat$treatment.assignment,
             prob = money.ri.probs,
             X = money.Xs)
money.Ys = genouts(Y = money.dat$allocated.diff,
                   Z= money.dat$treatment.assignment,
                   ate = 0)
##generate outcomes
distout <- gendist(money.Ys,money.ri.perms, prob=money.ri.probs, X = money.Xs)
money.Ys.ci = genouts(Y = money.dat$allocated.diff,
                      Z= money.dat$treatment.assignment,
                      ate = ate)
distout.ci <- gendist(money.Ys.ci,money.ri.perms, prob=money.ri.probs, X = money.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

##cohen's D for stadndardized effects
cohen = cohen.d(d = money.dat$allocated.diff,
                f = as.factor(money.dat$treatment.assignment))

##add outcomes to matrix
outmat['allocation','int.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 0])
outmat['allocation','seg.mean'] = mean(money.dat$allocated.diff[money.dat$treatment.assignment == 1])
outmat['allocation','ate'] = ate
outmat['allocation','p.value'] = disp$greater.p.value
outmat['allocation','cohens.d'] = cohen$estimate
outmat['allocation','N'] = nrow(money.dat)


#################################
##now do same thing for all other outcomes
#################################

#####attribution analysis
attribution.Xs = as.matrix(attribution.dat[,x.vars])
attribution.ri.perms = genperms(attribution.dat$treatment.assignment,
                                blockvar = attribution.dat$Block.ID)
attribution.ri.probs = genprob(attribution.ri.perms)
ate = estate(Y = attribution.dat$attribution.diff,
             Z= attribution.dat$treatment.assignment,
             prob = attribution.ri.probs,
             X = attribution.Xs)
attribution.Ys = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = 0)
distout <- gendist(attribution.Ys,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
attribution.Ys.ci = genouts(Y = attribution.dat$attribution.diff,
                         Z= attribution.dat$treatment.assignment,
                         ate = ate)
distout.ci <- gendist(attribution.Ys.ci,attribution.ri.perms, prob=attribution.ri.probs, X = attribution.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = attribution.dat$attribution.diff,
                f = as.factor(attribution.dat$treatment.assignment))

outmat['attribution','int.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 0])
outmat['attribution','seg.mean'] = mean(attribution.dat$attribution.diff[attribution.dat$treatment.assignment == 1])
outmat['attribution','ate'] = ate
outmat['attribution','p.value'] = disp$greater.p.value
outmat['attribution','cohens.d'] = cohen$estimate
outmat['attribution','N'] = nrow(attribution.dat)


####social analysis
diff.Xs = as.matrix(social.dat[,x.vars])
diff.ri.perms = genperms(social.dat$treatment.assignment,
                         blockvar = social.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = social.dat$social.diff.scale,
             Z= social.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = social.dat$social.diff.scale,
                  Z= social.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = social.dat$social.diff.scale,
                f = as.factor(social.dat$treatment.assignment))

outmat['social.perceptions','int.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 0])
outmat['social.perceptions','seg.mean'] = mean(social.dat$social.diff.scale[social.dat$treatment.assignment == 1])
outmat['social.perceptions','ate'] = ate
outmat['social.perceptions','p.value'] = disp$greater.p.value
outmat['social.perceptions','cohens.d'] = cohen$estimate
outmat['social.perceptions','N'] = nrow(social.dat)



###physical analysis
diff.Xs = as.matrix(physical.dat[,x.vars])
diff.ri.perms = genperms(physical.dat$treatment.assignment,
                         blockvar = physical.dat$Block.ID)
diff.ri.probs = genprob(diff.ri.perms)
ate = estate(Y = physical.dat$physical.diff.scale,
             Z= physical.dat$treatment.assignment,
             prob = diff.ri.probs,
             X = diff.Xs)
diff.Ys = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = 0)
distout <- gendist(diff.Ys,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
diff.Ys.ci = genouts(Y = physical.dat$physical.diff.scale,
                  Z= physical.dat$treatment.assignment,
                  ate = ate)
distout.ci <- gendist(diff.Ys.ci,
                   diff.ri.perms, 
                   prob=diff.ri.probs, X = diff.Xs)
disp =  dispdist(distout, ate = ate, display.plot = F)

cohen = cohen.d(d = physical.dat$physical.diff.scale,
                f = as.factor(physical.dat$treatment.assignment))


outmat['physical.perceptions','int.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 0])
outmat['physical.perceptions','seg.mean'] = mean(physical.dat$physical.diff.scale[physical.dat$treatment.assignment == 1])
outmat['physical.perceptions','ate'] = ate
outmat['physical.perceptions','p.value'] = disp$greater.p.value
outmat['physical.perceptions','cohens.d'] = cohen$estimate
outmat['physical.perceptions','N'] = nrow(physical.dat)


##format the matrix for output
outmat = as.data.frame(round(outmat,2))

outmat[,'cohens.d'] = outmat[,'cohens.d']*(-1)



##output Appendix Table S3
table.s11 = xtable(outmat,
                  digits = 3)
rownames(table.s11) = c('Allocation','Attribution', "Social Perceptions", "Physical Perceptions")
colnames(table.s11) = c('Integrated','Segregated','Beta','P','D','N')

print(table.s11,'TableS11.tex',
      type = 'latex',
      floating = F,
      booktabs = T)


cat('\n \n Experiment 2 Results, Non-Skeptical Subjects Only: \n')
print(outmat)

